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Phase field models have been applied in recent years to grain boundaries in single-component 
systems. The models are based on the minimization of a free energy functional, which is constructed 
j^jfj' phenomenologically rather than being derived from first principles. In single-component systems 

r-{ t the free energy is a functional of a "phase field", which is an order parameter often referred to as 

, the crystallinity in the context of grain boundaries, but with no precise definition as to what that 

term means physically. We present a derivation of the phase field model by Allen and Cahn from 
classical density functional theory first for crystal-liquid interfaces and then for grain boundaries. 
The derivation provides a clear physical interpretation of the phase field, and it sheds light on the 
parameters and the underlying approximations and limitations of the theory. We suggest how phase 
field models may be improved. 

o ' 

PACS numbers: 05.70.Np, 61.72.Mm, 71.15.Mb 

I. INTRODUCTION 

The classical density functional theory (DFT) advanced by Haymet and Oxtoby 1,2 was developed to describe 
crystal-liquid interfaces. In this paper we will show how this theory may be approximated in a transparent way to 
derive phase field models of crystal-liquid and grain boundary interfaces in single-component systems. We will also 
discuss how classical DFT may be applied to grain boundaries in single-component systems. At specified temperature 
and chemical potential, atoms at the boundary are free to rearrange themselves to accommodate a fixed change of 
crystal orientation on cither side of the boundary. In principle DFT provides an exact treatment of these defects at 
the atomic level within a grand canonical ensemble. The application of DFT to grain boundaries is of considerable 
interest in itself as an alternative to grand canonical Monte Carlo simulations and molecular dynamics simulations. 
Grain boundaries in single-component systems have been modeled using phase field approaches, most recently by 
■ Carter and coworkers 3,4 and Nestler and Wheeler 5 . All of these studies make use of an Allen-Cahn type free energy 
functional, obtained in the case of Carter et al. after first integrating out an "orientation field" . These models are 
based on the minimization of a free energy, which is constructed phenomenologically, referring to symmetries, relevant 
physical parameters, effective interactions etc., rather than being derived from first principles. The free energy is a 
00 . functional of a "phase field" , which is an order parameter often referred to as the "crystallinity" in the context of 
grain boundaries. The crystallinity is interpreted as a measure of the local degree of structural order, but with no 
C**** , attempt to define it more precisely. 

Our aim in this paper is to derive a phase field model for crystal-liquid interfaces and grain boundaries in single- 
component systems from classical density functional theory, providing not only a clear physical interpretation of the 
" i phase field, but also an exposition of the various approximations along the way. The mapping also enables us to 
identify limitations of phase field modeling of interfaces, and to suggest how it may be improved. 

We begin by outlining classical DFT and its application by Haymet and Oxtoby to crystal-liquid interfaces, noting 
key assumptions. While this derivation focuses on a thermodynamic self-consistency equation, applying the same 
procedure directly to the grand potential leads to the Allen-Cahn equation and subsequently to a phase field model 
for grain boundaries recently introduced by Carter et al.. 3 . An interesting by-product of the derivation of a phase field 
model for grain boundaries is the result that the energy and width of a boundary are both predicted to be reduced 
when there is a short common vector in the reciprocal lattices of both crystals. 
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II. CLASSICAL DENSITY FUNCTIONAL THEORY 



Following an integration over the momentum degrees of freedom, the canonical partition sum 6 in d dimensions for 
TV identical classical particles at positions rj.,. . . , r/v € fi, where fl is an arbitrarily large volume, interacting via a 
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potential Vn(ti, . . . , r n ) in an external potential U(r) is 



Z N = ^K- dN Jd d n...d d v N ex P ^-/3(y JV (r 1 ,...,r n ) + ^£7(r i ))^ , (1) 



where 



\2rmr) y ' 

is the dc Broglie thermal wave length. The Laplace transform of the canonical partition sum is the grand canonical 
partition sum 

oo 

Z = e^ N Z N . (3) 

This suggests the introduction of the dimensionless one-body local potential, not to be confused with the effective 
one-particle potential introduced later: 

u(v) = ^-PU(v). (4) 

The dimensionless grand potential, W = W[u] = — InZ, is then a functional of the local potential u only, apart from 
the implicit dependence on the interaction potential, and the temperature which is assumed to be constant. 

Introducing an operator /3jv(r; ri, . . . , r n ) = J2i ^( r — r i) reveals that the particle density, which is the ensemble 
average (•) of this operator, is given by functional differentiation of the grand potential with respect to the local 
potential 

p([u];r) = ( / 3 JV (r;ri,...,r n ))= -^^y- (5) 

This can be seen by rewriting f3pN — PY^i=i ^( r «) as Ja^r u(r)f>N(r; . . .) in the grand canonical partition sum. 
Taking a Legendre transform, the free energy 

T[p] = W[u]+ ( d d r' P {v')u(y') (6) 
Jn 

produces the conjugates of the cumulants of the density, in particular 

u{v)=u{[p);v) = -^-T[p\. (7) 

The free energy of an ideal, i.e. non- interacting system, where Vn = 0, can be integrated, 

F iA .[p]= [d d rp(r)(\n(A d p(r))-l) (8) 
Jn 

and immediately gives rise to the barometric formula p(r) id = h~~ d exp(u id . (r)). In an ideal system, the external 
potential necessary to produce a given density profile is found simply by inverting the barometric formula. For a 
given density profile p observed in an interacting particle system, one can calculate the local potential that would 
be needed in an ideal system, u id , to produce the same density profile. The difference between the local potential, 
u, operating in the interacting system, and the ideal local potential is the effective one particle potential C([p];r), 
defined as 

C([p];r) =u id .([p};v) -u([p];r) = ln(A d p(r)) -u([p];r). (9) 

In other words, an ideal system with local potential u + C has the same density profile p as the interacting system 
with local potential u. 

Similarly, one defines the excess free energy 



-*\p\ = r\p\-?M 



(10) 
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and finds by differentiation 
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Sp(r) 



$[p]=C([p];r). (11) 



Further differentiation of the effective one-particle potential produces higher order correlation functions. In particular, 
the direct correlation function ([/>]; r, r') is defined as 

c ' 2,( W ;r ' r ' )s 5»)*l^ < 12) 

Eq. (12) immediately suggests an expansion of the excess free energy about a reference state, which in the following 
will be the infinite homogeneous liquid with constant density po, 

4>[p] = $[ p = Po ] + fd d rC([p^p };r)(p(r)-p ) + l fd d rd d r> C^([p= p }; r, r')(p(r) - p )(p(r r ) - p ) + • . . (13) 
Jsi L Jsi 

We use this expansion below to derive features of the crystal phase from the reference system. 12 

Using the Ornstein-Zcrnikc equation 6 the Fourier transform c(k) of the two-point direct correlation function over 
a domain with volume V" in a translationally invariant system, see Eq. (23), can be directly related to the structure 
factor S(k) 

S(k) = (l-p V"c(k)y 1 , (14) 

so that the direct correlation function has a direct physical and experimental meaning. It is, in fact, the only link 
to experiment, which raises the question of how it is possible that a perturbation theory about a liquid can describe 
features of a crystal. For example, one might doubt that the tetrahedral structure of crystalline silicon could be 
predicted from the two point direct correlation function of liquid silicon, which has an average coordination number 
of 6. On the other hand, the functional ([p]; r, r'), i.e. the direct correlation function including its functional 
dependence on the density within the entire system, contains all information about higher order correlations through 
functional differentiation. In Eq. (13) is evaluated only at p(r) = po, so the full functional dependence is 

suppressed. A key assumption of the theory is that by a judicious choice of po the expansion in Eq. (13) can be made 
acceptable. This debate has appeared prominently in the literature 7 . On a more technical level, the functional Taylor 
series Eq. (13) may have a finite "radius of convergence" in p(r) — p , not extending beyond the liquid phase, not 
least because a phase transition introduces singularities. 



III. APPLICATION TO INTERFACES 

In all that follows the approximations are based on the two point direct correlation function C^ 2 \[p = po];r, r'), 
which is a function only of the separation |r — r' in a homogeneous system. To ease notation, the homogeneous, infinite 
reference liquid will be denoted by subscript 0, instead of explicitly carrying the argument p = p . In this notation, 
the one-particle potential of the reference system, which is assumed to be homogeneous and to have vanishing local 
potential u(r) = uq = j3p, is related to its density by 

Co = ln(A%) - Pp (15) 

consistent with Eq. (9). 

We use the reference system to calculate features of an interfacial system (distinguished by subscript i) for a given 
local potential M,(r), which amounts to finding a root pi(r) of m = 5Ti[p]/5p. Equivalently one can minimize the 
functional 

W([p], H) = HP] - [ A d r Ul {v)p{v) (16) 
Jsi 

with respect to p. At the minimum p = pi and W([/9i], [v,i]) reduces to the grand potential, as seen in Eq. (6). For a 
vanishing local potential in the interfacial system m = (3p. The expansion in Eq. (13) together with Eq. (15), Eq. (8) 
and Eq. (10), is the starting point for the DFT used by Haymet and Oxtoby 1 : 

W([p t ], [m = Pp]) = f dV [In (Mr')/p Q ) - 1] Pi(r') - <5 + / dV C oPo (17) 
Jsi Jsi 

-\[ dV / dV Cf(|r" - r'|) ( Pi (r') - p ) ( Pi (v") p ) . 
z Jsi Jsi 
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It is worth stressing that the only approximation made so far is the truncation of the functional Taylor series Eq. (13) 
at second order. 

To derive the results of [1,2] one requires that W differentiated with respect to pi vanishes. Using (r) = 

(2) (2) 

Cq (— r) = Cq (|r|) this produces the self-consistency equation 

= ln Pi (r)/p - /dVC< 2) (|r-r'|)( ft (r')-po), (18) 



n 



which can also be obtained by expanding the effective one particle potential in a functional Taylor series. It is very 
difficult to find a root pi for this equation satisfying particular boundary conditions far from the interface. 
Physically it is very appealing to represent the density in a pseudo-Fourier sum 1 ' 7 



p,(r)=p (l + ]>>n(rK k " r ] (19) 



with reciprocal lattice vectors k n indexed by n e Z d . For simplicity we have assumed the crystal has a monatomic 
basis. In crystals with more than one atom in the basis there will be additional phase factors in the pseudo-Fourier 
expansion. At a crystal-liquid interface there is only one crystal lattice, but at a grain boundary there is one set of 
reciprocal lattice vectors for each crystal, and the expansion in Eq. (19) is over the union of reciprocal lattice vectors 
of the two crystals. In the following derivation we will not make use of any specific properties of {k n } other than its 
discreteness and completeness (see Eq. (27)), as well as its orthogonality (see Eq. (24)), which all necessitate a finite 
Fourier domain, introduced in Eq. (20) as the volume V. For grain boundaries it follows that the derivation applies 
only to cases where a three-dimensional coincidence site lattice (CSL) 8 exists with a reasonably small unit cell. 

Owing to the spatial dependence of the coefficients jU n (r), Eq. (19) is not a true Fourier sum. However, it can 
accommodate every density Pi(r), even if the density is not periodic. While the spatial dependence has, a priori, the 
unfortunate consequence that individual coefficients /x n (r) cannot be calculated by projecting pi(r) onto the functions 
e -ik„r^ spatial dependence of the coefficients can be thought of as a means to separate length scales. The atomic- 
scale density variations are captured by the exponential exp(ik n r), while more gradual changes in density are captured 
by p a (r). We discuss below how the p n (j) for different r can be obtained by taking a local Fourier integral within a 
volume V(r), 

Pn(r) Po = V- 1 [ d d r Pl (r)e- lk " r for n^O (20) 
Jv(v) 

and po(l + / i o(?)) = V^ 1 J v ^d d r Pi(r). Here V(r) denotes the domain the integral is running over, a parallelepiped 
of volume V which is a unit cell of the coincidence site lattice or a multiple thereof. 

Using the new representation Eq. (19) in the self-consistency equation (18) and expanding the p n (r') about r gives 
for Eq. (18) 




v" Po J2< 



zk n r 1 T/"^_ \ A g£k n r 



p n {r)c(k n ) - z(V r • V k )M„(r)c(k n ) - i(V r • V k )V„(r)c(k n ) 



(21) 



dV^ 2) (|r-r'|)e- 4k(r - r,) =V"c(k) (22) 



n 



has been used, based on the definition 



c(k) = V"- 1 f d d r" e- 4kr "C< 2) (r") . (23) 
Jv" 

The central idea is that V" is so large that vanishes outside this volume and that Q is sufficiently large that the 
failure of Eq. (22) due to the shift by r when writing Eq. (22) as Eq. (23) is negligible. 

So far, only two approximations have been made, the functional Taylor series of the excess free energy and the Taylor 
series for the pseudo-Fourier coefficients /i n . It is perhaps not surprising that Eq. (21) remains just as intractable as 
Eq. (18), which becomes obvious by taking a Fourier transform over the domain V(r) centered at r G 51 on both sides 
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FIG. 1: Cartoon of a local Fourier transform. The full line shows the actual density profile p(r) within a system of size f2 
parameterized by an expression like Eq. (19) which contains a space dependent pseudo-Fourier coefficient. A Fourier transform 
taken over the small region V marked by the dashed lines around r, produces coefficients that reproduce the density locally 
(shown as thick line) and can be continued periodically throughout the system (dotted line). The inset shows the two point 
direct correlation function C' 2 ^(r) which drops to a negligible value outside the domain V" . 



of Eq. (21): 



V- 1 [ d d re- lk ™ r \n I 1 + V,! n (r)e' k ° r = 

Mr) V n / 

JV(t) „ 



M„(r)c(k n ) - i(V r • V k K(r)c(k n ) - - (V r • V k ) 2 /in(r)c(k n ) + . 



(24) 



The problem is the r-dependence of the coefficients /x n (r) which renders the orthogonality of the exponentials on the 
right hand side unexploitable. On the other hand, the representation in Eq. (19) has not yet been fully exploited, in 
particular, no use has been made so far of the many degrees of freedom in the parametrization. 

Haymet and Oxtoby achieved a separation of length scales by utilizing the parametrization Eq. (19). Firstly, they 
introduced a family of density profiles indexed by r, 



(25) 



which coincides with Eq. (19) for r = r. Secondly, they replace /z n (r) by /z n (r) in Eq. (21), requiring that the fi a (r) 
at fixed r are solutions of a slightly different problem, 



In l + ^ Mn (r)e 



ik„r 



,zk n r 



/i n (r)c(k n ) - ^(V r • V k K(r)c(k n ) - l(V r • V k ) 2 ^ n (r)c(k n ) + . . . 



(26) 
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and therefore 



V- 1 f d d re- 4kmr ln ( 1 + V Mn (F)e 4k " r ) 

J V(r) V n / 

V- 1 [ d d re-*™ r V"p Ve' 1 



'V(F) 



Mn(r)c(k n ) - i(V r • V k )M„(r)c(k n ) - l(V r • V k ) 2 M„(r)c(k n ) + . . . 
Mm(r)c(k m ) - i(V r ■ V k )/x m (r)c(k m ) - i(V r • V k ) 2 /i m (r)c(k m ) + . . . 



, (27) 



corresponding to Eq. (24). For a complete, discrete set of reciprocal lattice vectors {k n } Eq. (27) is equivalent to 
Eq. (26); in other words, Eq. (26) is implied by Eq. (27) only if (27) applies to the entire set of reciprocal lattice 
vectors. 

The solutions of Eq. (27) for each fixed r are periodic density profiles (see Fig. 1) of systems in a certain periodic 
external potential and it is plausible to assume that solutions in the form p n (r) exist. Moreover, Eq. (26) evaluated 
at r = r coincides with Eq. (21), so that any set of solutions of Eq. (27) also represents a solution of Eq. (24). In 
other words, Haymet and Oxtoby construct a family of physical systems parameterized by r in an external potential, 
the density profiles of which are given by p n (r) and can be used to construct a solution of the interfacial problem 
Eq. (24) by evaluating Eq. (25) at r = r. 

Because Eq. (27) can also be obtained by assuming in Eq. (24) the /z n (r) vary so little within V(r) they can be 
replaced by /x n (r) (see Eq. (20)), going from Eq. (24) to Eq. (27) is equivalent to a separation of length scales. This 
separation of length scales validates Eq. (20). It is illustrated in Fig. 1, where the Fourier coefficient /U n (f) obtained 
by taking a Fourier integral over a small domain centered at r, gives rise to a density profile that can be periodically 
continued throughout the system, but coincides almost perfectly with the actual density pi(r) within V(r). 

Small V are desirable so that the separation of length scales as encapsulated in Eq. (20) applies with no or only 
small corrections. For a grain boundary, the volume V is a unit cell of the CSL, so that large £ boundaries are 
expected to be less reliably treated. Similarly, small V" are desirable so that the expansion in Eq. (21) has only small 
corrections, which implies that shorter-ranged direct correlation functions are easier to handle in this formalism. The 
size of both volumes is to be compared to the scale on which /x n (r) change. 

We see in Eq. (27) that the separation of length scales embodied in Eq. (25) simplifies the original equation (24) 
dramatically by enabling the orthogonality of the functions exp(ik n r) to be exploited. Instead of finding roots fi a (r) 
for all k n and all r G 17 of the original integral equation Eq. (24) simultaneously, the integro-differential equation 
Eq. (27) can in principle be integrated, as both sides are local in r. As shown in ref. [2], the /U n (r) can be considered 
positions of particles in a (complicated) potential at "time" r projected on the interface normal. In DFT, the next 
steps would consist in determining the necessary degrees of freedom, incorporating all symmetries and prescribing a 
procedure to find the roots /x n (r) for every n and r. 

IV. DERIVATION OF PHASE FIELD MODELS 
A. Crystal-liquid interfaces 

We will now build on the approximations, expansions and parameterizations described in the previous section to 
derive an Allen-Cahn phase field model for crystal-liquid interfaces from the grand-potential Eq. (17). In the next 
sub-section we will derive a phase field model for grain boundaries. 

The grand potential VV in the form Eq. (17) is reparameterized using Eq. (19) and the integrals over Q rewritten 
^ /pdV = ^d-rV-^dV: 

W(\pi], [ui = = Jd d rV- 1 J v ^ d r > [\n ^ + E MOe*"'') - l^o ^1 + £ Mn(r'K k " r '^ - $o (28) 

+ f dV Co Po \pl [ d d r V-'f dV f dV C^\v" r') V fc(r> m (r")e' 
Jn z Jn Jvc?) Jn ~~z 



*(k n r'+k m r") 



As in the calculation of the Fourier coefficients of C^- 2 \ see Eq. (23), the equality between the single integral over Q 
and the double integral over individual volumes V(r) centered at r ignores surface terms and holds only in the limit of 
infinite or periodic SI, because the volume V(r) at a point r close to the surface of SI might not be fully within 0. After 
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expanding the coefficients /Um(r") about r' and using the definition Eq. (23) (in the notation Vfcc(— k m ) = Vfc|_ km c(k) 
and noting that c(— k) = c(k) etc. ) the last triple integral becomes 



V" fd d rV- 1 f dVY>„(r' 



Mm(r')c(k m ) - i(V r • V k )A* m (r')c(k m ) - -(V r • V k ) 2 Mm(r')c(k m ) + . . . 



a !(k n +k m )r' 



(29) 



which is very similar to the right-hand side of Eq. (21). Assuming a perfect separation of length scales, so that the 
/z n (r') do not change within the volume of a unit cell V, allows us to replace /z n (r') by /x n (r) within the integrals over 
V(r) and to make use of the orthogonality of the exponentials, so that W becomes 



= Jd^V^J d d r' (\n ^l + ^/i n (r)e ,k " r ' 



-1 ft 



i + ^>„(?K k " r ') 



(30) 



A*n(r)c(k n ) - i(V r • V k )^„(r)c(k n ) - -(V r • V k ) 2 M„(r)c(k n ) + . . . 



This is the final result for the grand potential. Differentiating with respect to /U-j(r) gives 
8 



<^-j(r) 



W(\pi],[ui = 0ii]) 



v- 1 
-pIv" 



f d d r' p e tk - jr ' In ( 1 + V p n (r)e 4k " r ' J 

J V{r) V n / 



(31) 



W(?)c(kj) - »(V r • V k )/Zj(f)c(kj) - -(V r • V k ) 2 W (f)c(kj) + . 



using, again, c(kj) = c(k_j), Vfcc(kj) = — Vfcc(k_j) etc. and k_j = — kj. Moreover, we made use of 

5 



Jdyf'(y)g(y) = -g'(x) 



(32) 



which is obtained by ignoring surface term in an integration by parts. This is particularly important when trying 
to include into a functional the first derivatives, V r /z n (r), reminiscent of a Rayleigh friction term in a mechanical 
interpretation. 

Eq. (31) is to be compared to Eq. (27), the original result by Haymet and Oxtoby, which is reproduced by Eq. (31) 
by requiring g ^ - ^ W = 0. The time evolution of /ij(r) is usually set equal to that derivative, 



-M 



W 



(33) 



with mobility M (which can be absorbed into the definition of time), driving the system to the above mentioned root. 
This is precisely the mechanism used in many phase field models. 

Further simplifications are needed to write the grand potential Eq. (30) in Allen-Calm form. The set of independent 
amplitudes jU n (f) is replaced by a single scaling amplitude (f>(r) by writing 



Mn(r) = 4>(r)vl 



(34) 



with constant amplitudes chosen to represent a crystalline phase, so that cf>(r) — 1 corresponds to the crystalline 
phase, while 4>(r) = suppresses all structure, corresponding to a liquid with density po- The field <f> therefore is 
to be interpreted as the order parameter, proportional to the local amplitude of the atomic density waves, i.e. (f> is 
the crystallinity. The common amplitude </>(r) is necessarily real, because p being real implies ^ n (r)* = p_ n (r) and 
together with = p°_ n (see below) we therefore have <f>(r)* = </>(?). 

Collecting all local contributions in the density w((j)(r)), the grand potential Eq. (30) simplifies to 



MM) 



d d rw((p{r)) 



(35) 



z(V r • V k )<K?Kc(k n ) + -(V r • V k ) 2 0(?)^c(k n ) 



8 



which can be simplified further by ignoring the surface terms of the integral J d d r (f>(r)\7 r (f>(r) so that the remaining 
integral containing derivatives of <j)(r) has the structure 

/ d d f 0(r) ]T M°_ n Mn(V r • V k ) 2 ^)c(k n ) • (36) 
Jn 

To simplify this expression, the set of reciprocal lattice vectors k n in the sum must be reduced. We use the index n° 
to indicate that we are considering only a subset {k n o} of all possible reciprocal lattice vectors {k n }. In the simplest 
approximation, the set is reduced to the set comprising only the shortest (non-vanishing) reciprocal lattice vectors, 
{k n o}. This set {k n o} forms a "star" which means that elements of this set are related by point symmetry operations 
of the reciprocal lattice (for further details see the Appendix, Subsection 1), so that the star is invariant under these 
operations. If the real-space lattice is FCC, for example, the reciprocal lattice is BCC and the eight nearest neighbor 
reciprocal lattice vectors 111, 111, 111 form the star of shortest length vectors. By including only the shortest 
reciprocal lattice vectors the symmetry and lattice constant of the crystalline structure are described correctly by 
the truncated Fourier expansion. The approximation can be systematically improved by including stars of longer 
reciprocal lattice vectors. 

Symmetry requires that all coefficients /z° as introduced in Eq. (34) associated with any of the vectors k n o within 
the same star have equal magnitude. Inversion symmetry in the real-space lattice ensures that H°_ n o = A^o, if we 
impose that the center of inversion coincides with a site. Since /J^ = ^°- n a by reality of p, all /i° are real and 
therefore equal, /i° = Mo- It is this equality which is needed to simplify Eq. (36), by allowing us to place n°_ n /z„ in 
front of the remaining sum, which is dealt with in the following. 

The direct correlation function of the bulk liquid is isotropic, c(k) = c(|k|), so that the second derivative with 
respect to wave-vector components 1 < a < d and 1 < (3 < d in Eq. (36) can be written as 

d ka d k0 c(V = ^c'(|k|) - ^5'(|k|) + ^c"(|k|) . (37) 

In a sum of the form 



^(V r • V k ) 2 0(?)c(k n o) = E E d ka d kfj d Sa d Sl Mr)c(Ko) 

k„o k„ a=l (3=1 

= EE fe^lM) - ^^(IM) + %^(IM)) d Sa d^(r) , (38) 

k n o a, n n I n I 

where the k n o run over all elements of the star {k n o}, the contributions of off-diagonal elements, a ^ (3, vanish after 
taking the summation over the star k n o, if one assumes a cubic crystal system. Only the diagonal elements remain 
and by symmetry each Cartesian component contributes J2{k i. fc 2 = (q/d)\k\ 2 , where q — |{k n o}| is the cardinality 
of the star, i.e. the number of elements in {k n o}. The dimension d is the dimension of the space spanned by the star. 
For example, the 8 nearest neighbors in a BCC lattice (the reciprocal lattice of an FCC lattice) would have q = 8 and 
d = 3. 

The proof of this simplification is based on the Great Orthogonality Theorem of group theory 9 as detailed in the 
Appendix. The assumption of cubic symmetry ensures that there is always a three-dimensional, unitary, irreducible 
representation of the point group. The assumption can be lifted if one allows for anisotropic terms (see the Appendix, 
Subsection 1), reflecting the anisotropy of the lattice. In the following, we consider only cubic crystal systems. 

The sum in the integrand of Eq. (36) now becomes 

^^^^(Vr-Vk) 2 ^)^^) = ^fcils' ( | ko |) + ^''(|k o |)) M ^(?)V 2 0(?) (39) 

where we have used |k | to denote the magnitude of any member of the star {k n o}, and we have introduced the 
coupling ek defined as 

£k » = - (nk^ 5 ' (|ko|) + ! 2 " (|ko|) ) ^ ■ (40) 



9 



The expression for ek simplifies further because the structure factor Eq. (14) usually peaks at the shortest reciprocal 
lattice vector 7 , so that c' vanishes and c" < 0. In that case 

e ko = ||c"(|k |)|>ug . (41) 

which enters the grand potential Eq. (35) as follows 

Wac(M) - Jfr (w(<f>(r)) + \p 2 V"e ko (V r 0(r)) 2 ) , (42) 

ignoring surface contributions again. Functional differentiation of Eq. (42) with respect to the phase field <j>(r) then 
reproduces the Allen-Cahn equation 10 

4> = -M ^V"||e"(|ko|)|/i§V^(r)) (43) 

in which the relaxational assumption 

(with mobility M) has been made to drive W[4>] to a minimum with respect to <f>. 

The above analysis may be improved by including more stars of reciprocal lattice vectors than just the shortest. 
Indeed, the entire reciprocal lattice can be decomposed into disjoint stars without double-counting. The sum over 
more than one star {k n o} produces an effective coupling 

e = E e Ko ( 44 ) 

{k„o} 

replacing 6k in Eq. (42). 



B. Grain boundaries 



As noted in Section III we confine ourselves to misorientations where there is a CSL with a relatively small three- 
dimensional unit cell to ensure that the finite Fourier domain V(r) in the separation of length scales (see Eq. (20)) 
does not lead to significant errors. In practice this limits the treatment to grain boundaries in cubic lattices, where 
such CSLs arise frequently 

In the simplest treatment of a grain boundary we consider one set of shortest length reciprocal lattice vectors in 
each crystal, which we call §; = {k^ } and § r = {k^ } for the left and right crystals respectively. Again, both these 
sets form stars, see Appendix, Subsection 1, and they are related by the rotation that generates the misorientation 
in the bicrystal. The atomic densities deep in the left and right crystals are described by density waves with wave 
vectors that are elements of these stars, and with equal amplitudes fjr n0 = fi and /x^ = Mo respectively, similar 
to the situation in the crystal-liquid interface. That does not mean, however, that all wave-vectors have non-zero 
amplitudes: For example, if a wave-vector in the right star is not also member of the left star, then its amplitude fi n0 
vanishes identically by symmetry. Amplitudes arc non-zero and equal within the respective stars. Again, it is this 
equality that is used in the following to simplify the expressions. 

In the spirit of Eq. (34), we approximate the spatial dependence of the amplitudes {/i n o(r)} by introducing a phase 
field 4>(r) for the grain boundary in the form 

Mn o (r) = (1 - </>(?)) mL« + <K?Ko • (45) 

To satisfy the boundary conditions far from the boundary plane <j)(r) must vary from deep in the left-hand crystal 
to 1 deep in the right-hand crystal. But this form of {^ n o(r)} does not allow for density waves which have non-zero 
amplitudes only in the grain boundary core, thereby restricting the degrees of freedom available to the system. There 
is no such restriction in place at the starting point of the derivation, Eq. (30), which included, in principle, all k- vectors 
of the direct sum of the two stars, i.e. the entire "DSC lattice" 8 . 

At closer inspection, the parameterization (45) of /i n o(r) makes a slightly different use of the "crystallinity" (j>(r), 
as it forces the system to be, on average, a superposition of both crystalline lattices: As one is gradually "switched 
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off", the other is gradually "switched on", so that </>(r) is the crystallinity of the right lattice and 1 — <j>(r) is the 
crystallinity of the left lattice. Reciprocal lattice vectors common to both lattices are predicted by Eq. (45) to be 
constant throughout the bicrystal, if /J = /Iq. Even if fi l ^ /1q the parameterization (45) disallows the possibility of 
a disordered grain boundary structure where all the {/z n o(r)} are locally zero. To resolve this problem, an alternative 
form of Eq. (45) is needed, for example [i n o = (1/2)0 ((1 + (p)^, + (1 — </>)Mn°) w ^ tn 4> varying from —1 to 1. Yet, 
by continuity such a form would force the grain boundary to be disordered somewhere 4 . In the following, we will use 
Eq. (45), which leads to the Allen-Cahn equation in a natural way. 

To see how the two different stars enter, the following derivation is presented in some detail. Inserting Eq. (45) in 
Eq. (30) produces the sum 



k 



[(i - m) mLo + <^K«] (v r • v k ) 2 [(i - m) mL" + <K?Ko] c(k° ) (46) 

which can be rewritten as four sums, 

Ek (W)/£o (Vr'Vk) 2 (W)/iLo 

+ E k " (W)/4> (v r -v k ) 2 4>^ 
+ E k „ <M r n « (v r -v k ) 2 (W)/4> 
+ E k „ <M r n « (v r • v k ) 2 w n0 . 

In the grand potential these sums are integrated over all space, see Eq. (36). A term of the form </>(V r • Vk) 2 (l — 4>) 
therefore equals (1 — </>)(V r • Vk) 2 following integration by parts and ignoring surface terms, which in turn equals 
— 0(V r • Vk) 2 after ignoring surface terms again. The first sum in (47), for example, then reads 

Ek n0 (W)/4> (v r -Vk) 2 (i-^Xo 

- Ek„o ^ (V r -Vk) 2 <f>n l n0 (48) 

= E§ ; <^L° ( V r • V k) 2 <ML° 

where the last equality is due to the coefficients /j l n0 vanishing for those k-vectors that are not member of the left 
star. For those k-vectors that are members, on the other hand, the corresponding coefficients /z^ are all equal to /j, l , 
which can be placed outside the sum. Applying the same procedure to all four sums in Eq. (47), gives 13 

E [(1 - 0(?)) mL" + 0(?Ko] (V r • V k ) 2 [(1 - </>(?)) n l n o + 0(?Ko] (49) 

s ; us r 

= ^ <^)( V r • V k)W) +MS 2 E ^( ? )( V r ' V O 2 0(?) - 2/i{,M5 E ^ ? K V r ' V k)V(?) . 

Si S r I 

The intersection I = §; nS r contains the reciprocal lattice vectors common to both stars, for which the product H^oH^o 
is non-zero. The intersection is assumed to be itself a star (see the Appendix, Subsection 1 for details). To Eq. (49) 
the simplification based on the Great Orthogonality Theorem (see the Appendix) can be applied. Because the space 
spanned by I is potentially only a sub- vector vector space of R d , a projection matrix Pi needs to be introduced, which 
projects any vector of R d to this sub-vector space. For the cubic systems we are focusing on here, this sub-vector 
space has either dimension dj = in which case Pi = 0, or dj = d in which case Pi = 1 is the identity matrix, or 
di = 1 in which case one spatial direction, say z, can be chosen to coincide with the common direction. Using the 
Great Orthogonality Theorem, Eq. (49) becomes 

^(r)) + 0(?Ko] (V r • V k ) 2 [(1 - </>(?)) + 0(?Ko] c(k„o) 



Eta 



S;US r 

= (jtf+tfW ffcil 5 '(|k |) + ^"(|k |) ! Vioiy) (50) 



- 2 M <^(?) ^2'(|k |)V 2 - ^L-5'(|k |)V r PlVr + |S"(|k |)VrPlVr j 0(F) , 

where the stars S; and § r each contain q elements, I contains qj elements, and d\ is the dimension of the space spanned 
by I. The two terms in Eq. (50) have the same prefactor if /J = fi^ = fi , which is a physically sensible choice we 
adopt henceforth. 
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Suppose there is no misoricntation between the crystal lattices. Then there is no grain boundary and Eq. (50) 
should reduce to zero. With no misorientation we would have S; = S r = I, so that qi = q, di = d and Pi = 1. 
Substitution of these values into the right hand side of Eq. (50) shows that it does indeed vanish. 

We identify in Eq. (50) an isotropic coupling term in which 

e> = ~ 2m ° (n^ s ' (|ko|) + ! 5 " (|ko|) ) + fe 2 ' (|kol) ) (51) 

multiplies — </>(r)V;: 4>(r). There is also an possibly anisotropic coupling term multiplying —(j)(r) V r Pi V r 4>(r): 

E ^ 2 " s (-|i« <i ' (|k ° l) + ^" (|k ° l) ) ' (52) 

This last term is not necessarily anisotropic, because for the particular choice of stars, Pi might be or the identity. 
In a cubic system, if Pi is non-zero and non-trivial for a given pair of stars, it is equal to any non-zero, non-trivial Pi 
produced by any other pair of stars. If the approximation is improved by adding further pairs of stars, they will all 
generate the same types of terms, namely either multiplying or V r PiV r . 
The resulting grand potential has again square gradient form 

Wac(M) - jfr (w(4>{t)) + \plV" (e; (V r 0(r)) 2 + e„V P #r)iWr0(r))) , (53) 

We note the difference in sign for the two couplings, Eq. (51) and Eq. (52), which becomes most obvious when c' = 0. 
Thus, we have shown that common reciprocal lattice vectors reduce the square gradient term in the interfacial energy, 
resulting in a lower interfacial energy and a smaller interfacial width, the characteristic scale of which is given by the 
coupling £j + e a , which has the dimension of a length squared. The smaller this coupling, the smaller the penalty for 
rapid changes in <f>. The amount of the reduction by common reciprocal lattice vectors depends on the magnitude of 
the second derivative c"(|ko|), which decreases with increasing |ko|. Thus the shorter the common reciprocal lattice 
vectors the larger their influence on the width and energy of the boundary. The influence of common reciprocal lattice 
vectors on the boundary energy has been known for a long time 8 , but we are not aware that their influence on the 
boundary width has been noted before. 

V. DISCUSSION 

In this paper we have shown how the classical density functional theory of Haymet and Oxtoby may be modified to 
produce an Allen-Cahn type free energy functional first for crystal-liquid interfaces and then for grain boundaries. For 
both types of interface the phase field is identified with the amplitudes of atomic density waves, providing a physical 
interpretation of the "crystallinity" of phenomcnological phase field models. 

Let us consider first the approximations that were required to derive an Allen-Cahn free energy functional from 
classical density functional theory. To obtain Eq. (21) only two approximations were made. Firstly, the excess free 
energy functional was expanded and the resulting series was truncated at the second order term, see Eq. (18). Secondly 
the Taylor series for £t n (r') about r in Eq. (21) was also truncated at the second order term. In principle both these 
expansions may be taken to higher order terms, although not without a significant increase in complexity. This would 
be the natural way to extend a phase field model to higher order terms with parameters related to fundamental 
quantities such as higher order correlation functions. 

The key step which lead to equation (30) for the grand potential W, and eventually to a phase field model, was the 
separation of length scales. However, this was not an approximation as we discussed after Eq. (27). 

The approximations mentioned so far were those of classical density functional theory of crystal-liquid interfaces 
as implemented by Haymet and Oxtoby. To obtain a phase field model for a crystal-liquid interface two further 
approximations were necessary. Firstly, the spatial dependence of all /x n (r) was assumed to be described by a single 
field, the "phase field" (f>(r). While classical density functional theory naturally handles a larger set of density waves 
with independent amplitudes phase field modeling in single-component systems to date has considered only a 
single order parameter, which we have identified as being proportional to the local amplitude of the atomic density 
waves. The obvious deficiency is the highly restricted nature of the configurational phase space made available to the 
crystal- liquid interface. When the same approach is adopted for a grain boundary, see Eq. (45), the restrictions are 
even more severe. Secondly, all /j l n0 and pL r nQ within a star are assumed to be identical and real, which suppresses any 
relative phase factor of the density waves. For example, if there were a rigid translation of one crystal relative to the 
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other on cither side of the grain boundary this would lead to complex density wave amplitudes in one crystal. The 
simplistic assumption made in Eq. (49) eliminates the possibility of any rigid body relative translation. 

It is clear that there are quite severe limitations of current phase field models of grain boundaries in single- 
component systems, at least when they are interpreted in the framework of classical density functional. Can we use 
DFT to indicate how these phase field models may be improved? We have already seen that going beyond the second 
order term in the expansion Eq. (18) will introduce higher order correlation functions and hence more information 
about structure and bonding. Indeed, going to at least third order correlation functions would seem to be necessary 
in covalent crystals such as silicon. This would introduce higher order terms in the grand potential Eq. (30) and in the 
Allcn-Cahn free energy functional Eq. (42) with parameters that are directly related to the higher order correlation 
functions. But perhaps the most obvious and pressing improvement would be the introduction of a second independent 
phase field to describe the amplitudes of two sets of density waves, one for each crystal. Far from the interface the 
crystal structures are related by a rotation. If we also wish to be able to predict a relative rigid body translation 
between the crystals, with components parallel and perpendicular to the interface, it will be necessary to allow the 
density wave amplitudes of one crystal to be complex. 
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In this Appendix, we derive the equations based on the Great Orthogonality Theorem, as used in Eq. (39) and (50). 
In particular, we will identify their prerequisites and possible extensions. 
The aim is to simplify expressions of the form 



where we may interpret M a p as an element of a matrix M. Introducing bra and ket notation for convenience the 
matrix M can be expressed as, 



In Eqns. (A.l) and (A. 2) the sums run over a set of vectors {k n o} c R d of equal length, which, for the time being, 
we will call a "star" . This term will be defined properly in the next section. 

The star is picked from the reciprocal lattice, and it is expected to display some point group symmetries of this 
lattice, in the sense that these point group operations effect only a permutation of members of the star, so that the 
star itself remains invariant. We assume that the star under consideration is invariant under a group of order <?, which 
has a unitary, irreducible representation {U n }, n = 1, . . . , g, and that this representation has the same dimension d as 
the vectors that make up the star. Acting with U n from the left and with t/t from the right on both sides of Eq. (A. 2) 
merely amounts to a permutation of the summands, because of the invariance of the star. Therefore M commutes 
with any group element, 
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APPENDIX: SUMMING OVER A STAR 




(A.l) 




(A.2) 



U n MUl = M 



(A.3) 



so that J2 n U n MUl = gM and therefore 
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d 




(A.4) 



where in the last equality we have used the Great Orthogonality Theorem (for unitary representations), 
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(A.5) 



n 
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By construction, the trace Tr M is the sum over the squares of the moduli of all vectors in the star, which all have 
the same length, which we call |k |. So, if the cardinality of the star is q, then TrM = <7|k | 2 and one finally arrives 
at the general result 



l/3 = ||k | 2 <W (A.6) 



provided there exists a group under which the star is invariant and for which there is a unitary, irreducible represen- 
tation of the same dimension as the vectors that make up the star. 

The fact that the Great Orthogonality Theorem applies only to irreducible representations of the group is a limitation 
of the above derivation, which motivates the following extension. By construction, the rank of the matrix M is d, 
which is the dimension of the vector space the k n o are taken from. A priori, it is unknown whether there exists a 
group under which the star is invariant and that has an irreducible, unitary representation with dimension d. Eq. (A.6) 
applies only if a suitable group and a suitable representation exists. 

If the star spans a sub vector space of dimension d' less than d, one cannot expect that it is invariant under 
a representation with dimension d. For example, the star Ki = {10 0,010,100,010} is not invariant under any 
irreducible three-dimensional representation of any group, because it will always contain some group elements that 
will rotate the star out of the xy-plane in which all members of the star are located. So, Eq. (A.6) cannot apply to 

On the other hand, considering only the two-dimensional sub vector space spanned by Ki (the xy plane), the 
star K2 = {10,01,10,01} derived from Ki by appropriate projection, is invariant under some unitary, irreducible 
representations of C\ v , some of which have dimension 2 (and some have dimension 1). Of course, the original star Ki 
is also invariant under some three-dimensional representations of C^ v , but none of them is irreducible. 

If a star is contained entirely in a vector space of dimension d' < d, its members can be expressed in rf-dimcnsional 
space as 

|k n o) = OB \k' n „) and (k n o| = (k^ | B+Ot (A.7) 

where O is an arbitrary rotation matrix of rank d, B is a d x d' projection matrix with Bij — Sij and |k^ ^ is a vector 

in M d . The matrix B increases the number of elements in the vector from d' to d and the rotation matrix O rotates 
the resulting vector to an arbitrary position. The vectors {|k^ )} also form a star, now in the sub- vector space M. d , 
so that 

£Ko}(k' n o| = |m| 2 l (A.8) 

if there exists a group under which the star {|k^ )} is invariant, and which has a (unitary) irreducible representation 
of dimension d! . Therefore 

^|k n o)(k n o| = ||ko| 2 (OB^Ot) (A . 9) 

using |k | = |k(,| and applying OB on Eq. (A.8) from the left and B^O^ from the right. The matrix (OBB^O^) = P 
(rotate, remove element, rotate back) projects a vector in R d onto the vector space isomorphic to R d . For example 

(A.10) 

for the star Ki introduced above, using its representation in the xy plane, K2. The construction of the vectors 
k n o from k^ in Eq. (A.7) represents an important limitation: At first sight, Eq. (A. 9) seems to apply to every 
star that has some degree of symmetry. However, not every star can be written in the form Eq. (A.7). For example, 
K3 = {101,011,101,011} has a d v symmetry, but this star spans the entire R d and so there is no star {|k^ )} C M d 
that would fulfill Eq. (A.7) with d' < d. 

In a more compact form, the general result is 




,afc„o,/3 = ||k | 2 P^ , (A.ll) 
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provided there exists an unitary, irreducible representation of a group under which the star {|k n o}} c K d is invariant, 
with a dimension d' identical to that of the space spanned by the star. Furthermore, if d! = d, then P reduces to the 
identity. Eq. (A. 11) can be used in the form 

EE fc »°.« fc »»^»^ = f l k o| 2 E P ^ 9 ^ = f M'V*^ (A.12) 

a/3 k n o a/3 

which leads to Eq. (50). 



1. Definition of a star 



So far in this Appendix a star has been defined as any set of A:- vectors of equal length. We now adopt the standard 
definition of a star, which will enforce the conditions required for Eq. (A. 6) or (A. 11) to apply: A star is generated 
from an initial vector by applying all symmetry operations in the form of a unitary, irreducible representation of the 
point group of the lattice in the appropriate coordinate system. 14 This representation is chosen so that the lattice 
itself is invariant under its action, so that the star represents a finite subset of the lattice. The star generated is highly 
degenerate, giving rise to the so-called "little groups" . 

There are two important caveats: Firstly, only cubic crystals have point group symmetries with three-dimensional 
irreducible representations. Secondly, the intersection of two stars does not necessarily have a irreducible with a 
dimension equal to that of the space spanned by the intersection. Generally, both points can be addressed in the 
same way, namely by decomposing stars or the intersection thereof into smaller stars which then have the required 
properties. But that produces additional projection matrices, so that, for example, the isotropic V~ terms in Eq. (39) 
and (50) split into multiple terms with different projections of the form VyPVy. 

All the cubic point groups, and only the cubic point groups, have three-dimensional, unitary, irreducible representa- 
tions, and the intersection of any two stars generated by one of these representations is either empty, one-dimensional, 
or it coincides with both stars. For cubic crystals all equations derived above apply without restrictions, producing 
an isotropic term and a anisotropic one with a single, preferred direction. 
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